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Numerical Study of Wave Propagation in Uniaxially Anisotropic Lorentzian Backward 
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The propagation and refraction of a cylindrical wave created by a line current through a slab of 
backward wave medium, also called left-handed medium, is numerically studied with FDTD. The 
slab is assumed to be uniaxially anisotropic. Several sets of constitutive parameters are considered 
and comparisons with theoretical results are made. Electric field distributions are studied inside and 
behind the slab. It is found that the shape of the wavefronts and the regions of real and complex 
wave vectors are in agreement with theoretical results. 

PACS numbers: 02.70.Bf, 41.20.Jb, 77.22.Ch, 77.84.Lf 



I. INTRODUCTION 

Metamaterials have received much attention during the 
last years, because they possess unusual electromagnetic 
properties, like, for example, the opposite directions of 
phase and group velocities. Double negative (DNG) ma- 
terials have negative permittivity and permeability, and 
they belong to the class of metamaterials. These me- 
dia that are capable of supporting backward waves, have 
been also called backward wave (BW) media in the lit- 
erature 1^. In BW media, the refraction phenomenon is 
anomalous in the sense that the power flow is refracted 
negatively, i.e. to the same side of the normal of the in- 
terface. As discussed in 1], it is not necessary for the 
medium to be a DNG medium to be able to support 
backward waves, because anomalous refraction can also 
be realized with anisotropic media with only one negative 
material parameter. We will use the term BW medium 
throughout the rest of the paper. 

The pioneering work on B W materials by Veselago 1^ , 
where slab lenses were mentioned, has gained much at- 
tention during recent years, despite some differing opin- 
ions on the subject 3]. Isotropic BW materials are often 
called Veselago materials. The possibility to realize a 
perfect lens with isotropic BW slabs was discussed by 
Pendry in 4]. Numerical and theoretical considerations 
of wave propagation in isotropic BW slabs excited with a 
line current above the slab were presented by Ziolkowski 
and Heyman in 0- Guidance of waves in a slab of uni- 
axially anisotropic metamaterial has been theoretically 
discussed by Lindell and Ilvonen in 0. 

Controversial opinions regarding negative refraction ef- 
fect and the perfect lens call for more detailed studies of 
the wave phenomena in backward-wave media. In this 
paper, wave propagation through uniaxially anisotropic 
BW slabs is numerically studied, and comparison is made 
with the theory presented in The theory predicts 
that there are regions in certain BW media, where the 



wave vector becomes complex, thus resulting in expo- 
nentially decaying waves. These regions are bounded by 
the asymptotes of the wave vector surfaces, which can 
be shown to be hyperbola. We study these phenom- 
ena numerically using the finite-difference time-domain 
(FDTD) method in a 2D-problem of a line current radi- 
ating in the vicinity of a BW slab. Also, the existence 
of surface waves on the interface between free space and 
BW medium is demonstrated with an example case. The 
BW medium is realized with Lorentzian constitutive pa- 
rameters having a single pole pair. 

An example problem with some theoretical discussion 
are presented in section 3. Results from the numerical 
simulations are shown and discussed in section 4. Our 
numerical simulations show that the wave propagation 
and refraction phenomena heavily depend on the param- 
eter choices of the BW medium and are qualitatively in 
agreement with the theory. 



II. NUMERICAL MODEL OF DISPERSIVE 
MEDIUM 

The constitutive relations for a frequency dispersive 
isotropic medium read 



D e(cj)E, B /i(cj)H. 



(1) 



Negative permittivity and permeability are realized with 
the Lorentz medium model. The expressions for the per- 
mittivity and permeability are of the form 
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(2) 



This model corresponds to a realization of BW materials 
as mixtures of conductive spirals or omega particles, as 
discussed in [3. In this artificial material both electric 
and magnetic polarizations are due to currents induced 
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on particles of only one shape, which provides a possibil- 
ity to realize the same dispersion rule for both material 
parameters, as in (|2j. Note that the medium realized 
by Smith et. al. is built using different principles 8]. 
For the uniaxial materials that we consider in this paper 
we assume that the negative components of the material 
parameters are realized by small uniaxial spiral inclu- 
sions (racemic arrays with equal number of right- and 
left-handed particles) and possess frequency dispersion 
defined by The positive components of the material 
parameters are equal to the free-space permittivity and 
permeability values, assuming that there are no particles 
oriented along these axes. 

Equations (QJ and Q form the basis of the used FDTD 
model for BW materials. The two most important known 
FDTD methods for modeling dispersive materials like 
Lorentz materials are the recursive convolution method 
and the auxiliary differential equation method, where the 
constitutive parameters are expressed with the help of 
susceptibility. In the first method, D and E and B and 
H are related through a convolution integral. This ap- 
proach is rather tedious. Another possibility is to use the 
auxiliary differential equation technique, which is slightly 
easier to implement. In this latter method, the polariza- 
tion current associated to each Lorentz pole pair is in- 
troduced. These two models are discussed in detail in 

i- 

A third method to discretize fields in Lorentz medium, 
classified as direct integration method in 10], is based 
on the direct discretization of the PDE representing 
the time-domain equivalent of the simplified frequency- 
domain constitutive relation. The proposed discretiza- 
tion scheme is a modification of this method. The idea 
is to transform (^) into the time domain using the rela- 
tion juo ^ d/dt with one integration before discretization 
using center differences. Usually, FDTD models based 
on the constitutive relation are directly (after multipli- 
cation with the denominator) discretized, as discussed 
in a summary of FDTD algorithms for dispersive media 
in To]. We found in 11] that one integration prior to 
discretization leads to much better accuracy and, what 
is also important, to considerably better stability proper- 
ties than the direct discretization without integration. In 
the numerical simulations, we have used the model which 
is discussed in detail in [llj . 



With these choices, we obtain e{(jj) = for all uj and 



-1 at a; 



The spatial reso- 



lution Ax = Ay = 1.5 cm is used throughout the sim- 
ulations. To be able to demonstrate the properties of 
BW materials, the incident spectrum is quite narrow, 
so that the relative constitutive parameters are close to 
minus one for the frequencies having significant spectral 
content. Absorbing boundary conditions are used to ter- 
minate the computational domain at the outer bound- 
aries of the lattice. For simplicity, we have used Liao's 
third order ABC, although more sophisticated ABC's are 
available. The use of usual ABC's requires a small gap 
between the outer boundary of the computational space 
and the BW material slab. The chosen coordinate axes 
and the problem geometry is shown in Figure^ Next, we 
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FIG. 1: The slab problem under consideration and the chosen 
coordinate system. 

briefly present some theoretical results from 1] that are 
important here for comparison purposes with our numeri- 
cal results. Wave propagation in BW slabs is studied with 
different value combinations of the medium parameters 
jj^x^ l-iy and e^. For our TE-polarized case, the Poynting 
vector Ste can be shown to read 1] 



TE 



2kor]ojiiojiixlJ^y 



\Eo\' 



(3) 



where \Eo\ is the amplitude of the TE-polarized electric 
field of the plane wave, k^^; is the wave vector with two 
cartesian components kx and ky^ ko = uo^cqiiq is the 

free space wave number, and r^o = y/W^ ^^^^ 
space wave impedance. Denoting the angle between the 
outwards-pointing unit normal vector \Xy and the wave 
vector \^TE by ^, one can derive the dispersion equation 
and solve it for the wave number as in JLJ. The result is 



III. AN EXAMPLE PROBLEM AND 
THEORETICAL DISCUSSION 



kTE{0) = ko. 



cos2 (9 + ^ sin^ 

My 



(4) 



Consider a z-directed line current in free space located 
at a distance dg from a BW slab of thickness d. Let the 
interface between free space and the BW slab be located 
at = 0. The problem space is two-dimensional, with 
the field components Hx , Hy , and Ez . The peak of the 
incident spectrum is at uOp = 5.0 • 10^ rad/s, and the pa- 
rameters in ((21) are the following: = ^Om = 1.0 • 10^ 
rad/s, LJ% = = 4.8 • lO^^ (rad/s)^, Fe = F^ = 0. 



For any given set of parameters, we may plot the pro- 
jection curves of the wave vector surfaces in x?/-plane. 
For certain choices of the medium parameters, the wave 
number becomes complex, resulting in exponentially de- 
caying waves inside the BW material. 

Any physical power flow must be directed downwards 
away from the source. This requires that 



• Ste < 0. 



(5) 
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The necessary condition for any transmission is 



^ <0. 



(6) 



Clearly, the phase velocity is directed oppositely to the 
power flow provided that /i^^ < 0. Negative refraction of 
the Poynting vector requires that 



fly < 0. 



(7) 



For the interface problem, we will also need to know 
conditions for the existence of surface waves. The input 
impedance on the surface filled by a uniaxial material is, 
for TE-polarized fields. 



TE 



(8) 



where the normal component of the propagation factor 
reads 



'TE 



^ {uj^e^fiy - kl) 



(9) 



with the square root branch defined by Ivci^Pte) < 0. kx 
is the propagation factor along the surface. Thus, the 
eigensolutions for an interface between this medium and 
free space satisfy the following equation: 



UJ/J.X 



DTE 







where 



(3o = Vc^^eo/io - k'^, 



(10) 



(11) 



with Im(/3o) < 0. 

Let us consider now conditions for the existence of sur- 
face waves along this interface. In this case both betas 
must be imaginary, of course with negative imaginary 
parts: /3te = —Jolte and (5^ = — jao, where aTE > 
and ao > 0. The eigenvalue equation becomes 



aTE 



0. 



(12) 



Obviously, if all the material parameters are positive, this 
equation has no solutions, but if fix < 0, surface wave 
solutions are possible. This is well known for interfaces 
with free-electron plasma. 

Five different cases will be considered in the follow- 
ing. In the first case, we choose jUx > 0, /i^ < 0, and 
Cz < 0. The wave vector surface is a two-sheeted hyper- 
boloid with axis parallel to the x-axis. The asymptotes of 
the hyperbola in x?/-plane divide the plane into regions of 
complex and real wave vectors. Waves that are propagat- 
ing parallel to y-ax.is are supposed to decay exponentially 
inside the slab. 

In the second case, /j^x > 0, /j^y < 0^ Cz > 0. The wave 
vector surface is a two-sheeted hyperboloid with the axis 



parallel to the surface normal. In our 2D-case, real wave 
vectors exist inside the region bounded by the asymptotes 
of the hyperbola associated to the wave vector surfaces. 

As a third case, we consider the situation complemen- 
tary to the second case in the sense that the signs of 
all the parameters are changed. Notice that this does 
not affect the shape of the wave vector curves. We have 
/jjx < 0, fly > 0^ and < 0. In fact, the Poynting vector 
is refracted positively in this case (see Q), but this is an 
interesting case anyway because of the aforementioned 
contrast with respect to the second case. 

The fourth case consists of the usual isotropic BW slab 
{fix = fly < 0^ and Cz < 0) where focusing and nega- 
tive refraction phenomena are present. The wave vector 
curves are ellipses (in our case of two equal parameters 
they are circles), and real wave vectors exist everywhere 
inside the slab. 

The fifth case is specially chosen to show the existence 
of surface waves in the case when fix < 0^ fiy < 0^ and 
Cz > 0. We can easily see from (0} that there are nei- 
ther real wave vectors nor backward waves. However, we 
have found that surface waves on the interface are easily 
excited in this case. Let us now present the numerical 
results. 



IV. NUMERICAL RESULTS AND 
COMPARISON WITH THE THEORY 



In the first four cases, we show the electric field dis- 
tribution at three suitable chosen increasing time steps 
to illustrate the wave propagation and refraction phe- 
nomena. In the fifth case, we illuminate a rectangular 
cylinder to see the surface waves. Whenever a constitu- 
tive parameter is said to be positive, it is supposed to 
be a constant and equal to the free space permittivity 
or permeability. Negative material parameters obey the 
Lorentzian dispersion rule and equal — eo or —fiQ at the 
center frequency. 



A. Case I: fix > 0, /i^ < 0, < 

In this case, the theory shows that the wave vectors 
are complex inside the slab within a region bounded by 
the asymptotes of a hyperbola. Inspection of Figure [2 c) 
reveals that there is indeed a region in the slab, where 
the electric field in negligible all the time. There are 
some fields within the slab just under the source. The 
hyperbola-shaped wavefronts propagate obliquely down- 
wards inside the slab and the power flow is refracted neg- 
atively. The distance of the source from the first interface 
is discretized with 10 cells, and the thickness of the slab 
corresponds to 80 cells. 
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(a) (b) (c) 

FIG. 2: (a) The electric field Ez penetrating into the slab. Slab boundaries are indicated by black lines, (b) The electric field 
within the slab. Notice the region in the center of the slab, where the field amplitudes are very small, (c) It appears that the 
wavefronts inside the slab are hyperbolas, in agreement with the theory. 



B. Case II: /i^ > 0, /i^ < 0, > 

In Figure [2l a) , a cylindrical wave is penetrating into 
the BW slab. In Figure 01 b), some numerical disper- 
sion is visible. There are significant fields inside the re- 
gion where the theory yields real wave vectors, while the 



fields are rather small elsewhere. Despite some dispersive 
effects, the wavefronts of constant field value are reminis- 
cent of hyperbolas. The phase velocity is directed down- 
wards. Some weak focusing of the power flow is seen in 
FigureOlb). The phase velocity inside the slab is directed 
downwards, as can be seen from the theory. 




(a) (b) (c) 

FIG. 3: (a) The electric field Ez penetrating into the slab, (b) The electric field within the slab. Some waves have already 
passed through the slab. It is seen that there are small fields in the region, where the wave vector is complex, as predicted by 
the theory, (c) Dispersive effects are more clearly seen inside the slab. The electric field is concentrated to the lower side of the 
slab. The wavefronts behind the slab outside the region of large amplitudes are prolate ellipses, ds = lOA^/, d = 80 Ay. 



C. Case III: /j^x < 0, /iy > 0, < 

To complete the analysis, we change the signs of the 
parameters of the second case. Clearly, the wave vec- 
tor surfaces as defined by are not changed. In fact. 



the Poynting vector is refracted positively in this case. 
However, the phase velocity is directed upwards [see 
and Q]. This phenomenon is clearly seen during the 
simulation. From Figure ^ we see that the electric field 
distributions inside the slab are quite similar to those of 
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Figure El except that the wavefronts are less distorted in are seen to be oblate ellipses with the center on the lower 
Figure^ In Figure ^c), the wavefronts behind the slab interface of the slab. 




FIG. 4: (a) The electric field Ez penetrating into the slab. Wavefronts are dramatically bent, (b) The electric field within 
the slab. Some waves have already passed through the slab. The wavefronts are seen to be hyperbolas, (c) Interestingly, the 
wavefronts behind the slab appear to be ellipses with the center on the lower interface of the slab, ds = lOA^/, d = 80 Ay 



D. Case IV: an Isotropic Slab with fix < 0, fiy < 0, 

ez <0 

Here we consider the usual isotropic BW (or double 
negative) slab with all the relative material parameters 
close to minus one. This case has been studied, for 
example, by Ziolkowski and Heyman in in the case 
of Drude slabs. We obtained quite similar results with 
our alternative discretization technique in the case of 



isotropic Lorentz medium. The electric field distributions 
are shown in Figure O We can calculate the positions of 
the foci from the slab thickness d and the distance dg of 
the source from the interface. Notice that we must have 
ds < dto have a focus inside the slab. The foci should ap- 
pear at 7/ = —ds inside the slab and at = —{2d — ds) be- 
hind the slab. The appropriate derivations can be found 
in |5|]. 




FIG. 5: (a) The electric field Ez penetrating into an isotropic BW slab, (b) The electric field within the isotropic BW slab. 
Some waves have already passed through the slab. The first focus inside the slab is seen in the figure, (c) The second focus 
behind the BW slab has become visible, ds — 20Ay, d — 60Ay 
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From Figure [51 a) we see that the electric field is con- 
centrated in the expected position of the first focus. In 
Figure Elc), the second focus behind the slab is also visi- 
ble. It takes some time for the foci to develop, as is seen 
from Figure inib), where the second focus is not yet seen. 
The incident spectrum has some small components for 
which the relative material parameters are not exactly 
minus unity. Hence the wavefronts are not perfect cir- 
cles as predicted by the theory. Anyway, these results 
are in agreement with the theory concerning negative re- 
fraction of the Poynting vector. However, no "perfect" 
focusing has been observed, meaning that the focus area 
is always not smaller that about half wavelength. Steady 
state solutions for the foci were not obtained. This same 
observation was also made by Ziolkowski and Heyman in 

i. 

E. Case V: /ia, < 0, /i^ < 0, > 

For this set of material parameters, the theory predicts 
that the waves decay exponentially everywhere inside the 
slab. However, we have found that surface waves are eas- 
ily excited in this case, in accordance with the theoretical 
prediction, see surface. To see this phenomenon, we il- 
luminate a rectangular cylinder with a pulse having a 
slightly broader spectrum. The electric field distribution 
induced on the surface of the cylinder is shown in Fig- 
ure El Figure El a) shows how the surface waves begin 



FIG. 6: (a) The surface waves begin to develop. The source 
is still on. (b) The source has been switched off, and surface 
waves exist around the cylinder. As is known from the theory, 
there are negligible fields inside the cylinder. 

to develop. The source is still on. In Figure El b), the 
source above the cylinder has been switched off. The 
surface waves have propagated to the opposite side of 
the cylinder as well. In agreement with the theory, there 
are no fields inside the cylinder except for the immediate 
vicinity of the surface of the cylinder. 



V. CONCLUSIONS 

Wave propagation and refraction phenomena in uni- 
axially anisotropic BW slabs have been studied with the 
FDTD method. Special attention was paid to the shape 
of the wavefronts and to the regions inside the slabs, 
where the wave vector becomes complex, thus resulting 
in exponentially decaying waves. The numerical results 
for anisotropic BW slabs were seen to qualitatively agree 
well with the theoretical results. The effects of nega- 
tive refraction, (imperfect) focusing, and surface wave 
excitation have been demonstrated. Potentially useful 
transformations of wave fronts between spherical, ellip- 
tical, and hyperbolical can be realized in homogeneous 
uniaxial backward wave slabs. 
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